Nonlinear analysis of the shearing instability in granular gases 



in 



R. Soto, M. Mareschal 
CECAM, ENS-Lyon, 46 Allee d'ltalie, 69007 Lyon, France 
M. Malek Mansour 

O ; Universite Libre de Bruges, Campus Plaine, Bvd. du Triomphe, CP 231, 1050 Bruges, Belgium 

o 

£NJ ■ It is known that a finite-size homogeneous granular fluid develops an hydrodynamic-like insta- 

| bility when dissipation crosses a threshold value. This instability is analyzed in terms of modified 

. hydrodynamic equations: first, a source term is added to the energy equation which accounts for 

■ the energy dissipation at collisions and the phenomenological Fourier law is generalized according 

to previous results. Second, a rescaled time formalism is introduced that maps the homogeneous 
cooling state into a nonequilibrium steady state. A nonlinear stability analysis of the resulting equa- 
tions is done which predicts the appearance of flow patterns. A stable modulation of density and 
temperature is produced that does not lead to clustering. Also a global decrease of the temperature 
' is obtained, giving rise to a decrease of the collision frequency and dissipation rate. Good agreement 

^ ' with molecular dynamics simulations of inelastic hard disks is found for low dissipation. 

a' 

> 

' I. INTRODUCTION 

a ' 

The understanding of the dynamics of granular fluids is crucial for various industrial processes. This has led to 
many investigations where theory, experiments and simulations are used in order to construct a predictive theory 
. There is hope that for moderate densities and slightly dissipative grains, grain dynamics may be described on 
large scales by using fluid hydrodynamics with slight modifications. 
£j ' One way to guess which changes to make in standard hydrodynamics is to use the tools developed by statistical 
mechanics in order to derive hydrodynamic equations. This is justified by the fact that a simple grain model, the 
inelastic hard sphere model, has been shown to reproduce most of the phenomena occuring in granular systems: in 
some sense, it has proven to contain the essential ingredients necessary to predict the peculiar physics observed |)] [uj. 

The scheme used in the study of nonequilibrium fluids can then be extended to granular fluids: "microscopic" 
simulations permit to compute the equation of state and values for transport coefficients which are then fed into the 
i guessed macroscopic equations. Comparison between direct nonequilibrium simulations of microscopic and macro- 
scopic models allows then to test the validity of the proposed macroscopic equations. This approach was used recently 
by two of us to investigate heat transport (heat being identified with the kinetic energy associated with the grains' 
motion) and it has been shown that Fourier's law has to be generalized with a density gradient term appearing in the 
expression for the heat flux Jll| . 

In the Inelastic Hard Sphere model (IHS), grains are spherical hard particles with only translational degrees of 
freedom. The energy dissipation is included through a restitution coefficient r lower than one. As for hard spheres, 
I . the collision is an instantaneous event and the grain velocities after a collision are given by 

§; vi=vi +-(l + r)[n-(v 2 -vi)]n (1) 

O ■ 1 

v 2 =v a --(l + r)[fi-(v a -vi)]n (2) 
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where n is the unit vector pointing from the center of particle 1 towards the center of particle 2. It is convenient to 
define the dissipation coefficient q = (1 — r)/2 which vanishes when collisions are elastic. In what follows units are 
chosen such that the disk diameter a and the particles masses are set to one. 

The IHS model, like the elastic hard sphere model, does not have an intrinsic energy scale. This means that two 
systems with the same initial configuration and with the particle velocities of one system being equal to those of the 
other multiplied by one scaling factor will follow the same trajectory but at different speeds. This lack of an intrinsic 
energy scale also implies a simple temperature dependence for all hydrodynamic quantities which may be found by 
simple dimensional analysis. 

When a fluidized granular medium is let to evolve freely in a box with periodic boundary conditions, the energy 
decreases continuously in time and the system remains homogeneous. This non-steady state is called the Homogeneous 
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Cooling State (HCS) and it is the reference state from which perturbations are studied: it is analogous to the 
thermodynamic equilibrium state for elastic systems Jl(].|l(i|,|l7| . In the IHS model this state is particularly simple and 
the energy decreases obeying the Haff 's law Q 

E(t) = (3) 
(l + t/t ) 2 

In order to avoid the cooling of the system, the simulations are made at constant energy: at every collision between 
the grains, the dissipated energy is redistributed by scaling all the velocities. This is equivalent to a time rescaling 
and can be treated as such in the equations for the continuous model. Using appropriate hydrodynamic equations, 
this homogeneous state is predicted to be unstable under certain conditions of density, system size and dissipativity 
P, ^0|JlS|Jl^ , |20| , pT|] . Considering the dissipativity coefficient q as a bifurcation parameter, while increasing q at constant 
density and number of grains, the system first develops an instability characterized by two counterflows, the shearing 
instability, and then either a clustering regime in which the density becomes inhomogeneous or a vortex state where 
many small vortices develop throughout the system. 

For given values of total number of grains N and number density p = N/V the shearing instability appears when 
the dissipation coefficient is larger than a critical value. In the low density and large system limit,the critical value is 
given (in two dimensions) by p2| 

Note that in the thermodynamic limit the system is always unstable for any finite dissipation. 

In this article we will study the shearing instability using a nonlinear hydrodynamic approach. In Sec. || we 
will develop a formalism that allows us to treat the HC S a s a non equilibrium steady state, and we will present a 



stability analysis around the steady state. Next, in Sec. Ill we will present the nonlinear analysis of the instability, 



obtaining expressions for the hydrodynamic fields beyond the instability. It will also be shown that the presence of 



the instability modifies the collision rate and energy dissipation rate. Finally, in Sec. IV we compare the predictions 
of the continuous model with molecular dynamics simulations of inelastic hard disks. Conclusions are presented in 
Sec. 0. 

In what follows we will treat the two dimensional case, but the extension to three dimensions is direct. 



II. RESCALED TIME FORMALISM 



Consider a two dimensional system composed of N grains interacting with the IHS model in a square box of size 
L that has periodic boundary conditions in both directions. Units are chosen such that Boltmann's constant and the 
particle masses are set to one. Granular temperature is defined analogously to the kinetic definition for classical fluids 



^£^-v) 2 (5) 



T 

N 

where v is the hydrodynamic velocity. 

The shearing instability has been predicted by linear analysis of the hydrodynamic equations for the IHS. Also an 
analysis of the first stages of the nonlinear regime has been done in Ref. |23j ]. As we shall see, the relative simplicity 
of the model allows for a complete nonlinear analysis as well. 

The hydrodynamic equations for the low-density IHS are similar to the usual hydrodynamic equations for fluids 
except that there is an energy sink term and the heat flux has contributions from the density gradient 22 24|,[l6|,[l5|] . 
The equations read 

^ + V • (pv) = (6) 

P (^ + (V ' V)V ) =- V]P W 
*° + ( v ' V ) T ) =-V-J--P: Vv-w (8) 
with the following constitutive equations 
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P l3 = P T5 l3 - Vo VT (|| + || - (V • v)5«) (9) 

J = -koVTVT - fj, Vp 

P 

uj = lu qP 2 T 3 / 2 

where 770, fc , po, and ojq do not depend on density or temperature but on the dissipation coefficient q. In particular, 
ujq and /iQ vanish with q. 

The total energy dissipation rate can be computed from the hydrodynamic equations, obtaining 

wdV (11) 



In the HCS, where density and temperature are homogeneous and there is no velocity field, the energy dissipation 
rate is 

^ = ^V^plT^ (12) 

This continuous cooling down has the disadvantage that any perturbation analysis must be done with respect to 
this non-steady state. To overcome this difficulty we propose a differential time rescaling ds = 7 dt such that in 
this rescaled time, energy is conserved. The transformation corresponds to a continuous rescaling of all the particle 
velocities such that the kinetic energy remains constant. This rescaling does not, however, introduce new phenomena 
since as we have mentioned the IHS does not have an intrinsic time scale, and thus a time rescaling gives rise to the 
same phenomena viewed at different speed. 

In order to keep the rescaled energy constant the appropriate value of 7 is given by 



The rescaled time hydrodynamic fields transform as 

v = 7V (14) 

T = 7 2 T (15) 

Pij = l 2 Wj (16) 

J=7 3 J (17) 

lo = 7% (18) 

where the tilde denotes the rescaled variable. 

Due to the non-constant character of the time rescaling, extra source terms appear in the hydrodynamic equations, 
which can be simplified using the relation 



1 d-y 1 



7 ds 2E(0) 
Suppressing the tildes everywhere, the equations now read: 



>dV (19) 



^+V-(pv) = (20) 
p(S + (W)v)=-V.iP + ^/^ 

-°(fj + (v " V)T ) =-V-J-iP:Vv- W +^|a,dF 
Note that the constitutive relations remain unchanged under time rescaling. 
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In the rescaled time, the HCS reduces to a non-equilibrium steady state with a continuous energy supply that 
compensates the energy dissipation. As there is no energy scale, we fix the reference temperature to be one, so that 
E(0) = N. The HCS is then characterized by {p = pH, T = 1, v = 0}. To study the stability of this state, we first 
introduced the change of variables: 

p = p H + Sp (21) 
v = Sv (22) 
T = 1 + 5T (23) 

Taking the (discrete) Fourier transform of the linearized hydrodynamic equations around HCS, it is easy to check 
that the transverse velocity perturbation decouples from the rest and satisfies the equation 

PH^^- = Xk6v k± (24) 



ds 



with 



AfcEE-^^o + ^p ; {k x ,ky} = 0,±1,±2, ... (25) 

For small values of luq (proportional to the dissipation) A& < 0, so that the perturbations 5~Vk ± decay exponentially 
(the case k — should not be taken into account since the center of mass velocity is strictly zero) . But there exists 
a critical value of luq for which A& vanishes, thus indicating the stability limit of the corresponding mode. The first 
modes to become unstable correspond to |k| = 1 (i.c, k x = ±1, k y = and k x = 0, k y = ±1). The instability 
threshold for luq is then given by 

Wo = -pW m 

The stability for the other modes has been studied previously fl6|| . In this last reference, it was shown that for, low 
dissipation, the first instability that arises is indeed the transverse velocity instability. 

The origin of the instability can be understood also in terms of the real time hydrodynamics. In fact, the relaxation 
of the transverse hydrodynamic velocity is basically due to the viscous diffusion, which depends on the system size, 
whereas the cooling process is governed by local dissipative collisions. There exists therefore a system length beyond 
which the dissipation of thermal energy is faster that the relaxation of the transverse hydrodynamic velocity. The 
latter will then increase, when observed on the scale of thermal motion, thus producing the shearing pattern. 

III. NONLINEAR ANALYSIS OF THE INSTABILITY 

In this section we propose to work out the explicit form of the velocity field beyond the instability. The calculations 
are tedious and quite lengthy, so that, here, we only describe explicitely the basic steps. We start by taking the 
Fourier transform of the full nonlinear hydrodynamic equations, obtaining a set of coupled nonlinear equations for 
the modes { 5pk, 8Tk, Svk^ , 6vk ± }, where Svk^ represents the longitudinal component of the velocity field. Close to 
the instability (ujq lDq, |k| = 1), the modes 5vi ± exhibit a critical slowing down since Ai ~ 0, whereas the other 



hydrodynamic modes decay exponentially (cfr. ref. |16|). On this slow time scale, i.e. s « 0(A|f *), the "fast" modes 
{Spk', 5T},] Sv^ ; 5vfc ± , k ^ 1} can be considered as stationary, their time dependence arising mainly through 8vi ± . 
Setting the time derivatives of these fast modes to zero one can express them in terms of the slow modes 5v± ± . If 
now one inserts the so-obtained expressions into the evolution equation for the slow modes, one gets a set of closed 
nonlinear equations for (adiabatic elimination p5| , p6| ), usually referred to as normal form or amplitude equations 
[G. Nicolis, Introduction to Nonlinear Science, Cambridge University Press (1995) ]. Note that such a calculation is 
only possible close to the instability threshold, where the amplitude of 5vi ± approaches zero as w — > ujq. On the 
other hand, in this limit the only Fourier modes that will have large amplitudes are the transverse velocity modes 
with wavevector equal to ± 2tt/L (i.e. |k| = 1). There are four such modes: 
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^= , C_ ni'r (27) 
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After some tedious algebra, one finds 



dAi 

^ H ds 
dA 2 

PH—r- 

ds 

where 



X1A1 - 4ir 2 L- 2 V0 (Ci \A^ 2 + C 2 |A 2 | 2 ) Ai (28) 
AiA 2 - 47r 2 i- 2 r? (Ci \A 2 \ 2 + C 2 |A| 2 ) A 2 (29) 



8(fc - no) - 3t? 
8(fc - Mo) - 2ryo 

2(fc -/i ) +7?p , , 

2(fc - Mo) - % 

The amplitude equations ( |28| ) and (|2^) admit three different stationary solutions: 

(a) Ax = 0, A 2 = (32) 
(6) 1^1=1, A 2 =0 (33) 
A x = 0, \A 2 \ = A (34) 



"» '^C^G 2 (35) 

where we have set 

Note that the phases of the above stationary solutions are arbitrary (recall that Ai and A_i are complex conjugate). 

The trivial solution (a) corresponds to a motionless fluid, whereas the solutions (b) are shearing states with the 
corresponding fluxes oriented either in the y or x direction. The mixed mode solution (c) represents a vortex state 
with two counter-rotating vortices in the box. 

Below the critical point (Ai < 0), the trivial solution (a) is the only stable one. As we cross the critical point 
this solution becomes unstable. A linear stability analysis of the Eqs. |2^ shows that the shearing states are stable 
provided 

£ > 1 (37) 

while the mixed mode solution is unstable. Satisfying Eq.|37] depends on the values of the transport coefficients. For 
small dissipativity it is always fulfilled, provided the number density remains relatively low. In fact, a mixed mode 
state has been observed recently in a highly dense system [ flO| |. 

The occurence of either shearing states depends on the initial state. In a statistical sense, they are equally probable. 
For example, in the case that the system chooses the solution {A\ — A, A 2 = 0}, the velocity field reads 

v = 21cos(2ttx/ j L) y (38) 

where y represents the unit vector in Y direction. The temperature and density perturbations are then given by 

ST = -{A) 2 [I + (1 — Ci) cos(47ra;/L)] (39) 
Sp = p ff (l) 2 (f - Ci)cos(47T2;/L) (40) 

The instability of the transverse velocity field thus gives rise to modifications of the temperature and density fields. 
The temperature decreases globally, since part of the energy is taken by the convective motion. Moreover, because of 
the viscous heating, the temperature profile exhibits a spatial modulation, i.e. it is higher where the viscous heating 
is higher. The density profile also shows a spatial modulation that keeps the pressure homogeneous (recall that for 
a two dimensional low density gas, dp « Ph^T + <5p, in system units). This modulation of the density, that was first 
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observed by Goldhirsh and Zanetti 0] (Fig. 2 of their article) in the study of the clustering instability, is stable and 
does not lead to clustering. 

Using the above expressions for the hydrodynamic fields, the energy density profile reads 




Assuming the molecular chaos hypothesis, the mean collision rate, and dissipation rate, u>, can be computed as 
well: 



These relations show that the global decrease of the temperature leads to corresponding decreases of the collision 
frequency and dissipation rate. 

It is important to note that the origin of the nonlinear coupling of slow modes lies in the viscous heating term and 
the state-dependence of the transport coefficients, and not in the usual convective derivatives ((v- V)v). In fact, as we 
have shown above, the shearing state produces a variation in the density and temperature fields that modify locally 
the value of the transport coefficients. This effect, which is negligible in classical fluids, can become very important in 
granular fluids, mainly because of the the lack of scale separation of the kinetic and hydrodynamic regimes. Contrary 
to normal fluids, here, the convective energy is comparable to the thermal energy. In fact, as will be shown in the MD 
simulations, in a well developed shearing state up to half of the total kinetic energy corresponds to the convective 
motion. The microscopic source of this phenomenon lies in the fact that only the relative energy is dissipated in 
binary collisions, i.e. the center of mass energy is conserved. In other words, the thermal energy is dissipated but the 
convective one is conserved. This asymmetry in the energy dissipation mechanism is at the very origin of the shearing 
instability. 



For the molecular dynamic simulations we have considered a system made of N — 10000 hard disks with a global 
number density pn = 0.005. Inelastic collision rules are adopted, with a dissipativity varying from o = 0.0 to o = 0.12. 
The boundary conditions are periodic in all directions. A spatially homogeneous initial condition is adopted, with 
velocities sorted from an equilibrium (zero mean velocity) Maxwellian distribution. We note that the density is low 
enough so that the system remains within the low-density regime. 

The simulations have been performed in the rescaled time. Computationally, this is achieved by doing a normal 
IHS simulation (event driven molecular dynamics [p7| ) , but at each collision the value of the kinetic energy is updated 
according to the energy dissipated. The instantaneous value of 7 is computed from the kinetic energy, allowing to 
evaluate all the rescaled quantities. Note that the s-time can be integrated in the simulation because 7 is a piecewise 
constant function, thus allowing to make periodic measurements in the system. Finally, to avoid roundoff errors, a 
real velocity rescaling is performed whenever the kinetic energy decreases by a given amount (typically 10~ 7 of the 
initial value). 

In each simulation, the collision frequency and temperature dissipation rate are computed with respect to the 
s-time, after the system has reached a stationary regime. In Fig. yj, the collision frequency and the temperature 
dissipation rate are presented as a function of the dissipation coefficient q. As expected, the shearing instability is 
associated with an abrupt decrease of these functions. It must be noted that the decreasing of the collision frequency 
is more that 30%, which corresponds to a global decrease of the temperature of more than 50%. This means that 
when the shearing is fully developed, about half the total kinetic energy is taken by the macroscopic motion. This 
phenomenon is typical of granular media and has no counterpart in classical fluids. 

To measure the critical point, qo, we fit the collision frequency according to the following piecewise function 




(42) 



(43) 



IV. MOLECULAR DYNAMICS SIMULATIONS 




(44) 



obtaining 



G 



q = 0.0686 (45) 

a = 0.0178 (46) 

ay = -0.167 (47) 

a 2 = 1.03 (48) 



Similarly, the dissipation rate is fitted according to 



M q < 9o 

b Q q + bi{q- q )q + b 2 (q - q ) 2 q q > qo 



(49) 



Using q a = 0.0686, one finds 



b Q = 0.000166 (50) 
by = -0.00466 (51) 
b 2 = 0.0546 (52) 

To compare these results with our theoretical predictions we need the explicit form of the transport coefficients, up 
to critical dissipativity. Unfortunately, there are no known expressions for them in the case the 2d IHS model in the 
low-density regime. However, as the critical dissipativity is small we can use the the quasielastic approximation for 
the transport coefficients (that is, taking the first non-trivial order in q) pSlpa 



uj = ^\fnq (53) 

k = (55) 
Mo = (56) 



The critical dissipativity and the amplitude of the shearing state are given by 




n= 7^ (57) 
A = PHh {W q (58) 

(59) 

For the presented simulation, the predicted critical dissipativity is 

g = 0.0628 (60) 

which shows a discrepancy of 8% with the observed value. This difference is consistent with the adopted approxima- 
tions. 

The predicted values for clq, ay, b , by (cfr. Eqs. E3, |57], and ^9|) are 

a Q = 0.0177 (61) 

ay = -0.146 (62) 

b = 0.000177 (63) 

by = -0.00438 (64) 

which are also consistent with the adopted approximations. 

Since the system is periodic, the developed convective pattern can diffuse in the direction perpendicular to the 
flow (the phases of the complex amplitudes are arbitrary due to Galilean invariance). As a result, the average 
hydrodynamic fields remain vanishingly small, mainly because of "destructive" interference. To overcome this difficulty 
we have performed another series of simulations, keeping periodic boundary conditions in the vertical direction, while 
introducing a pair of stress-free and perfectly insulating parallel walls in the horizontal direction (in a collision with 
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a wall the tangential velocity is conserved whereas the normal one is inverted). As a consequence, the total vertical 
momentum is conserved, which will be simply set to zero initially. 

The nonlinear analysis for this case is similar to the periodic one, except that, here, the direction of the flow pattern 
remains always parrallel to the walls. Furthermore, the unstable wavevector is now k = ir/L, because of to the fixed 
boundary conditions. As a result all the previous predictions remain valid, except that everywhere L must be replaced 
by 2L. 

We have used the very same number of particles and density for this series of simulations, but, of course, the different 
boundary conditions produce a new critical dissipativity. Performing the same analysis as before, the measured critical 
point and fit parameters turn out to be 



while the predicted ones are 



3o 
a 
at 
a 2 
bo 
h 
b-2 



(I 

Ol 



0.0163 
0.0178 
-0.562 
8.93 

0.000175 
-0.0179 
1.06 



0.0157 

0.0177 

-0.583 

0.000177 

-0.0175 



(65) 
(66) 
(67) 
(68) 
(69) 
(70) 
(71) 



(72) 
(73) 
(74) 
(75) 
(76) 



Eqs. M gg, and once L is replaced by 2L, indicate that the perturbation of the transverse momentum density 
(j = pvjhas a wavevector equal to k x — tt/L, while the density and energy density have wavevectors k x = 2n/L. In 
the simulations we computed the amplitudes of these Fourier modes using the microscopic definitions for the particle, 
momentum, and energy densities. Figs. |[ |^, and ^ show the predicted and computed Fourier modes amplitudes. The 
predictions are in good agreement with the simulations in the neighborhood of the critical point, showing that not 
only the average quantities like the collision frequency are well predicted, but also the whole hydrodynamic picture is 
correct. 



V. CONCLUSIONS 



Taking advantage of the lack of energy scale in the IHS model, a rescaled time formalism was introduced that allows 
us to study the homogeneous cooling state as a nonequilibrium steady state. Using a hydrodynamic description for 
granular media written with a rescaled time variable, the shearing instability has been studied in the nonlinear regime. 
It has been shown that the shearing state is the stable solution and its amplitude has been computed. The appearance 
of the velocity field produces that part of the kinetic energy goes from the kinetic to the hydrodynamic scale. In usual 
fluids this redistribution of the energy is negligible, but in granular fluids it can represent an important fraction of 
the total energy. This phenomenon is a manifestation of a global property of granular fluids: there is not in general a 
clear separation between the kinetic and the hydrodynamic regimes. This could lead to put into question the validity 
of a hydrodynamic description. Nevertheless, at small values of the dissipativity coefficient, predictions based on the 
nonlinear hydrodynamic equations are in excellent agreement with molecular dynamics simulations. Both the value of 
the critical dissipativity and the behavior after the instability has developed are well predicted. This is a remarkable 
result which shows again how robust the hydrodynamic fluid equations are when they are tested at time and length 
scales where their validity could be questioned. 
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FIG. 2. Amplitude of the k x — 2n/L cosine component of the density field as a function of the dissipativity q. The dots 
are the results of molecular dynamics simulations, the dashed line is the full nonlinear theory prediction, and the solid line is 
the nonlinear theory prediction where the critical point is taken from the fit. The simulations were done in a channel with 
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FIG. 3. Amplitude of the k x — 2n/L cosine component of the energy density field as a function of the dissipativity q. Graph 
symbols and simulation parameters are explained in Fig. 0. 
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FIG. 4. Amplitude of the k x = ty/L cosine component of the momentum density field as a function of the dissipativity q. 
Graph symbols and simulation parameters are explained in Fig. ti. 
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